#******************************************************************************
#******************************************************************************

# import libraries

# standard

import math

# external

import numpy as np

from scipy.optimize import fsolve

# local

from .classes import Pipe, InsulatedPipe

from ..common.fluids import Fluid, FluidDatabase

from ..common import formulas as core

# TODO: allow users to select whether to respect correlation validity limits

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

# constants

# laminar flow regime

FLOW_REGIME_LAMINAR = 'laminar'

# turbulent flow regime

FLOW_REGIME_TURBULENT = 'turbulent'

# transition flow regime

FLOW_REGIME_TRANSITION = 'transition'

# minimum reynolds number

LAMINAR_LOWER_RE_LIMIT = 1

# highest reynolds number where laminar flow rules still apply

LAMINAR_UPPER_RE_LIMIT = 2300 

# lowest reynolds number where turbulent flow rules still apply

TURBULENT_LOWER_RE_LIMIT = 3000

# highest reynolds number where turbulent flow rules still apply

TURBULENT_UPPER_RE_LIMIT = 1e8

# isothermal surface condition

CONDITION_ISOTHERMAL_SURFACE = 'isothermalsurface'

# constant heat flow condition

CONDITION_CONSTANT_HEAT_FLOW = 'constantheatflow'

# Colebrook White equation for the Darcy friction factor

DFF_COLEBROOK_WHITE = 'Colebrook-White'

# Haaland equation for the Darcy friction factor

DFF_HAALAND = 'Haaland'

# (First) Pethukov equation for smooth tubes

DFF_PETHUKOV = 'Pethukov'

# Nikuradse equation

DFF_NIKURADSE = 'Nikuradse'

# Blasius equation

DFF_BLASIUS = 'Blasius'

# Nikuradse-Prandtl-von Karman  equation

DFF_PRANDTL_NIKURADSE_KARMAN = 'P-N-K'

# Gnielinski nusselt number correlation

NU_GNIELINSKI = 'Gnielinski'

# Colburn nusselt number correlation (smooth tubes)

NU_COLBURN = 'Colburn'

# Dittus-Boelter nusselt number correlation (smooth tubes)

NU_DITTUS_BOELTER = 'Dittus-Boelter'

NU_MODELS = [
    NU_GNIELINSKI,
    NU_COLBURN,
    NU_DITTUS_BOELTER]

# dict with reynolds number limits for nusselt number correlations

NU_RE_LIMITS = {
    NU_GNIELINSKI: (3e3, 5e6),
    NU_COLBURN: (1e4, math.inf),
    NU_DITTUS_BOELTER: (1e4, math.inf),
    }

# dict with prandtl number limits for nusselt number correlations

NU_PRANDTL_LIMITS = {
    NU_GNIELINSKI: (0.5, 2e3),
    NU_COLBURN: (0.7, 160.0),
    NU_DITTUS_BOELTER: (0.7, 160.0),
    }

# dict with reynolds number limits for darcy friction factor models

DFF_RE_LIMITS = {
    DFF_COLEBROOK_WHITE: 
        (TURBULENT_LOWER_RE_LIMIT, TURBULENT_UPPER_RE_LIMIT),
    DFF_HAALAND:
        (TURBULENT_LOWER_RE_LIMIT, TURBULENT_UPPER_RE_LIMIT),
    DFF_PETHUKOV: (1E4,1E6), # Source: Cengel, 2003, page 441
    DFF_NIKURADSE: (TURBULENT_LOWER_RE_LIMIT, TURBULENT_UPPER_RE_LIMIT),
    DFF_BLASIUS: (2.32e3, 1E5), # (3e3, 2e5) according to  Avic et al. (2019)
    DFF_PRANDTL_NIKURADSE_KARMAN: (TURBULENT_LOWER_RE_LIMIT, 1E5) #(2e3,1e5)
        }
 
# Prandtl and Blasius: page 180
# Böswirth L, Bschorer S, Buck T. Technische Strömungslehre: Lehr- und Übungsbuch

# dict with relative pipe roughness limits for darcy friction factor models

DFF_RPR_LIMITS = {
    DFF_COLEBROOK_WHITE: (0,0.1),
    DFF_HAALAND: (0,0.1),
    DFF_PETHUKOV: (0,0), # smooth tubes only
    DFF_NIKURADSE: (1E-6,0.05),
    DFF_BLASIUS: (1E-6,0.05),
    DFF_PRANDTL_NIKURADSE_KARMAN: (1E-6,0.05)
        }

# list of all models for the Darcy friction factor

DFF_MODELS = [DFF_COLEBROOK_WHITE, 
              DFF_HAALAND, 
              DFF_PETHUKOV, 
              DFF_NIKURADSE,
              DFF_BLASIUS,
              DFF_PRANDTL_NIKURADSE_KARMAN]

# list of explicit methods

DFF_EXPLICIT_MODELS = [DFF_HAALAND, 
                       DFF_PETHUKOV, 
                       DFF_NIKURADSE,
                       DFF_BLASIUS]

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def FlowRegime(reynolds_number,
               re_laminar_lower=LAMINAR_LOWER_RE_LIMIT,
               re_laminar_upper=LAMINAR_UPPER_RE_LIMIT,
               re_turbulent_lower=TURBULENT_LOWER_RE_LIMIT):
    
    if reynolds_number < re_laminar_lower:
        
        return math.nan
    
    elif (reynolds_number <= re_laminar_upper):
        
        return FLOW_REGIME_LAMINAR
    
    elif (reynolds_number >= re_turbulent_lower):
        
        # turbulent region
    
        return FLOW_REGIME_TURBULENT
    
    else:
        
        # transition region
        
        return FLOW_REGIME_TRANSITION
    
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def MinimumReynoldsNumberCompleteTurbulence(friction_factor):
    
    return (
        (200/3.7+2.51)*
        (10**(1/(2*np.sqrt(friction_factor))))/
        np.sqrt(friction_factor)
        )
            
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorExplicitBlasius(reynolds_number):
    
    return 0.3164*(reynolds_number**(-0.25))

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorBlasiusEquation(reynolds_number, friction_factor):
    
    return (friction_factor - 0.3164*(reynolds_number**(-0.25)))
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

# Nikuradse = PrandtlNikuradseKarman

def DarcyFrictionFactorExplicitNikuradse(pipe_internal_diameter,
                                         pipe_absolute_roughness):
    
    # for rough pipes
    
    # source: Avic et al. (2019)
    # https://doi.org/10.1016/j.euromechflu.2019.07.007
    
    return (
        1.74-
        np.log10(4)-
        2*np.log10(pipe_absolute_roughness/pipe_internal_diameter))**-2
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorNikuradseEquation(friction_factor, d_int, e_eff):
    
    return (
        friction_factor-(1.74-np.log10(4)-2*np.log10(e_eff/d_int))**-2
        )
            
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorExplicitPethukov(reynolds_number):
    
    # source:
    # secondary:  Cengel 2003, p. 441, eq. 8-65
    
    # for smooth pipes
    
    return (0.790*np.log(reynolds_number)-1.64)**-2
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorPethukovEquation(reynolds_number, friction_factor):
    
    return (
        friction_factor-(0.790*np.log(reynolds_number)-1.64)**-2
        )

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorExplicitHaaland(reynolds_number,
                                       absolute_roughness,
                                       pipe_diameter):
                                   
    return (
        (1.8*np.log10(
            ((absolute_roughness/pipe_diameter)/3.7)**1.11+
            6.9/reynolds_number))**-2
        )

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorHaalandEquation(reynolds_number,
                                       friction_factor,
                                       absolute_roughness,
                                       pipe_diameter):
                                   
    return (
        1/np.sqrt(friction_factor)+
        1.8*np.log10(
            ((absolute_roughness/pipe_diameter)/3.7)**1.11+
            6.9/reynolds_number)
        )

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorPrandtlNikuradseKarmanEquation(reynolds_number,
                                                      friction_factor):
    
    return (
        1/np.sqrt(friction_factor)
        -
        2*np.log10(reynolds_number*np.sqrt(friction_factor))
        +
        0.8
        )

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorColebrookWhiteEquation(reynolds_number,
                                              friction_factor,
                                              absolute_roughness,
                                              pipe_diameter):
    
    return (
        1/np.sqrt(friction_factor)+
        2*np.log10(
            (absolute_roughness/pipe_diameter)/3.7+
            2.51/(reynolds_number*np.sqrt(friction_factor))
            )
        )

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactorEquationTurbulentFlow(
        reynolds_number,
        friction_factor,
        absolute_roughness,
        pipe_diameter,
        model=DFF_COLEBROOK_WHITE):
    
    if model == DFF_COLEBROOK_WHITE:
        
        return DarcyFrictionFactorColebrookWhiteEquation(
            reynolds_number, 
            friction_factor, 
            absolute_roughness, 
            pipe_diameter)
    
    elif model == DFF_HAALAND:
        
        return DarcyFrictionFactorHaalandEquation(
            reynolds_number, 
            friction_factor, 
            absolute_roughness, 
            pipe_diameter)
    
    elif model == DFF_PETHUKOV:
        
        return DarcyFrictionFactorPethukovEquation(
            reynolds_number,
            friction_factor)
    
    elif model == DFF_NIKURADSE:
        
        return DarcyFrictionFactorNikuradseEquation(
            friction_factor,
            pipe_diameter,
            absolute_roughness)
    
    elif model == DFF_BLASIUS:
        
        return DarcyFrictionFactorBlasiusEquation(reynolds_number,
                                                  friction_factor)
    
    elif model == DFF_PRANDTL_NIKURADSE_KARMAN:
        
        return DarcyFrictionFactorPrandtlNikuradseKarmanEquation(
            reynolds_number,
            friction_factor)
    
    else:
        
        # unknown method
        
        print('This friction factor model '+model+' is not known.')

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def DarcyFrictionFactor(reynolds_number,
                        pipe: Pipe,
                        flow_regime = None,
                        model = DFF_COLEBROOK_WHITE,
                        model_is = DFF_HAALAND,
                        tentative_solution = None,
                        ignore_model_input_ranges = True):
    """Calculates the Darcy friction factor for a given pipe and flow."""
    
    if flow_regime == None:
        
        flow_regime = FlowRegime(reynolds_number)
    
    if flow_regime == FLOW_REGIME_LAMINAR:
        
        # laminar flow

        return 64/reynolds_number
    
    elif (flow_regime == FLOW_REGIME_TURBULENT or
          flow_regime == FLOW_REGIME_TRANSITION):
        
        # turbulent flow
            
        # check if an initial solution has been provided
        
        if tentative_solution == None:
            
            # no initial solution has been provided, therefore:
            # a) use explicit methods
            # b) use explicit methods to generate an initial solution
            
            if model_is in DFF_EXPLICIT_MODELS:
                
                # get a solution using an explicit model
                
                _, tentative_solution = DarcyFrictionFactorExplicit(
                    pipe,
                    model_is,
                    flow_regime,
                    reynolds_number)
                
                # return tentative solution if the models match
                    
                if model_is == model:
                    
                    return tentative_solution
                
            else:
                
                print('The initial solution model is not explicit.')
                
                return np.nan
            
        #**********************************************************************
        #**********************************************************************
           
        # either the function has returned or tentative_solution is available,
        # either due to it having been supplied or generated
        
        # declare wrapper function
        
        def friction_factor_equation_wrapper(f):
            
            # if f < 0:
                
            #     return np.nan
            
            # else:
        
            return DarcyFrictionFactorEquationTurbulentFlow(
                reynolds_number,
                f,
                pipe.e_eff,
                pipe.d_int,
                model=model)
        
        # solve the system
        
        f =  fsolve(friction_factor_equation_wrapper, tentative_solution)
        
        # check the solution
        
        if f[0] > 0:
            
            # solution appears to make sense
            
            return f[0]
        
        else:
            
            # the solution does not make sense
            
            print('no friction factor could be produced')
            
            return np.nan
            
        #**********************************************************************
        
    else:
        
        # undefined flow regime
    
        print('unknown flow regime, no friction factor could be produced')
        
        return np.nan
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def InsulatedPipeSpecificResistance(h_internal,
                                    h_external,
                                    extra_specific_resistance,
                                    pipe: InsulatedPipe):
    
    if h_internal != 0 and h_external != 0:
    
        return (1/(math.pi*pipe.d_int*h_internal)
                +
                math.log(pipe.d_ext/pipe.d_int)/
                (2*math.pi*pipe.k)
                +
                math.log(pipe.d_ins/pipe.d_ext)/
                (2*math.pi*pipe.k_ins)
                +
                math.log(pipe.d_cas/pipe.d_ins)/
                (2*math.pi*pipe.k_cas)
                +
                1/(math.pi*pipe.d_cas*h_external)
                +
                extra_specific_resistance)
    
    else:
        
        # there is no heat transfer between them and the resistance is infinite
        
        return math.inf

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def PipeHeatTransferCalculations(pipe: Pipe,
                                 fluid_in: Fluid, 
                                 fluid_bulk: Fluid, 
                                 fluid_speed,
                                 heat_transfer_condition,
                                 temperature_surface=None,
                                 thermal_resistance=None,
                                 heat_flux=None):
     
    #**************************************************************************
    #**************************************************************************
    
    temperature_inlet = fluid_in.temperature # kelvins
    
    volumetric_flow_rate = core.VolumetricFlowRateFromMeanFluidSpeed(
        fluid_speed, 
        pipe.InternalCrossSectionalArea())
    
    mass_flow_rate = core.MassFlowRateFromVolumetricFlowRate(
        volumetric_flow_rate, 
        fluid_bulk.mass_density)
    
    # check which condition is valid
    
    if heat_transfer_condition == CONDITION_ISOTHERMAL_SURFACE:
        
        pipe_efficiency = (
            1-((temperature_inlet-temperature_surface)/temperature_inlet)*
            (1-math.exp(-1/
                      (fluid_in.specific_heat*
                       mass_flow_rate*
                       thermal_resistance)))
            )
        
        temperature_outlet = (
            temperature_surface-(temperature_surface-temperature_inlet)*
            np.exp(-1/
                   (fluid_in.specific_heat*
                    mass_flow_rate*
                    thermal_resistance))
            )
        
        # # lmtd
        
        # lmtd = LogMeanTemperatureDifference(temperature_inlet,
        #                                     temperature_outlet,
        #                                     temperature_surface)
        
        # heat transfer rate
        
        heat_transfer_rate = core.HeatTransferRateMassFlow(
            temperature_inlet,
            temperature_outlet,
            mass_flow_rate,
            fluid_bulk.specific_heat)
        
    else: # constant heat flux condition
    
        # heat flow rate
        
        heat_transfer_rate = heat_flux*pipe.InternalSurfaceArea()
    
        # outlet temperature
        
        temperature_outlet = (
            temperature_inlet+
            heat_transfer_rate/(mass_flow_rate*fluid_bulk.specific_heat)
            )
        
        # efficiency
        
        pipe_efficiency = (
            1+
            heat_transfer_rate/(
                mass_flow_rate*
                fluid_bulk.specific_heat*
                temperature_inlet)
            )
    
    #**************************************************************************
    #**************************************************************************
    
    return pipe_efficiency, temperature_outlet, heat_transfer_rate
                      
    #**************************************************************************
    #**************************************************************************

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def SpecificPressureLossInPipe(fluid_mass_density,
                               pipe_internal_diameter,
                               fluid_speed,
                               friction_factor):
    """Function for the specific pressure loss for a given flow rate and pipe."""
    
    return (friction_factor*
            fluid_mass_density*(0.5)*
            (fluid_speed**2)/pipe_internal_diameter)

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def HeatTransferRateInPipe(fluid_speed,
                           fluid_in: Fluid,
                           fluid_out: Fluid,
                           pipe: Pipe,
                           fluid_bulk: Fluid = None):
    """Function for the heat flow rate for twin pipes at a mean fluid speed."""
    
    if fluid_bulk == None:
    
        return (
            (fluid_out.specific_heat*
             fluid_out.mass_density*
             fluid_out.temperature
             -
             fluid_in.specific_heat*
             fluid_in.mass_density*
             fluid_in.temperature)
            *
            ((pipe.d_int/2)**2)*np.pi*fluid_speed
            )
        
    else:
        
        mass_flow_rate = (
            np.pi*((pipe.d_int/2)**2)*fluid_speed*fluid_bulk.mass_density)
    
        return core.HeatTransferRateMassFlow(fluid_in.temperature, 
                                             fluid_out.temperature, 
                                             mass_flow_rate, 
                                             fluid_bulk.specific_heat)

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def HeatFlowRateViaLMTD(heat_transfer_coefficient,
                        log_mean_temperature_difference,
                        surface_area):
    
    return (heat_transfer_coefficient*surface_area*
            log_mean_temperature_difference)

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def HydrodynamicEntranceRegionLength(reynolds_number,
                                     pipe: Pipe,
                                     flow_regime):
        
    # Source: Cengel2003, page 425
    # Primary sources:
    # Kays and Crawford (1993), Ref. 13
    # Shah and Bhatti (1987), Ref. 25
    # Bhatti and Shah (1987), Ref. 1
    # Zhi-qing (1982), Ref. 31

    if flow_regime == FLOW_REGIME_LAMINAR:
        
        # laminar flow
        
        # Source: 
        
        # W. M. Kays and M. E. Crawford. Convective Heat and
        # Mass Transfer. 3rd ed. New York: McGraw-Hill, 1993.
        
        # R. K. Shah and M. S. Bhatti. “Laminar Convective
        # Heat Transfer in Ducts.” In Handbook of Single-Phase
        # Convective Heat Transfer, ed. S. Kakaç, R. K. Shah, and
        # W. Aung. New York: Wiley Interscience, 1987.
        
        return 0.05*reynolds_number*pipe.d_int
    
    elif flow_regime == FLOW_REGIME_TURBULENT:
        
        # turbulent flow
            
        # M. S. Bhatti and R. K. Shah. “Turbulent and Transition
        # Flow Convective Heat Transfer in Ducts.” In Handbook
        # of Single-Phase Convective Heat Transfer, ed. S. Kakaç,
        # R. K. Shah, and W. Aung. New York: Wiley Interscience, 1987.
        
        #W. Zhi-qing. “Study on Correction Coefficients of
        # Laminar and Turbulent Entrance Region Effects in
        # Round Pipes.” Applied Mathematical Mechanics 3 (1982), p. 433.
        
        return 1.359*(reynolds_number**0.25)
    
    else:
        
        return math.nan

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def MechanicalPumpingPower(volumetric_flow_rate, pressure_drop):
    
    return volumetric_flow_rate*pressure_drop

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

# get a valid solution (f,Re) with the least amount of information possible
 
def DarcyFrictionFactorExplicit(pipe: Pipe,
                                model=DFF_HAALAND,
                                flow_regime=FLOW_REGIME_TURBULENT,
                                reynolds_number=TURBULENT_UPPER_RE_LIMIT,
                                ignore_model_input_ranges=True):
    
    # TODO: define exceptions for the hyhetra package
    
    # Returns valid (friction factor, Reynolds number) tuples using explicit
    # methods. This works for both laminar and turbulent flows but is arguably
    # more relevant for the latter, namely to generate initial solutions.
    
    # Laminar flows: reynolds number will be used if compatible
    # Turbulent flows: smooth pipe, only requires reynolds number
    # Turbulent flows: rough pipe, requires pipe information
    # Turbulent flows: any pipe, requires reynolds number and pipe information
    
    # model must be explicit, except if the flow is laminar
    
    if ((model not in DFF_EXPLICIT_MODELS) and 
        (flow_regime is not FLOW_REGIME_LAMINAR)):
        
        # # error: unknown explicit model
        
        # Warning('The Reynolds number does not match the flow regime.')
        
        # print('The Reynolds number does not match the flow regime.')
        
        return np.nan, np.nan   
    
    #**************************************************************************
    
    # if flow regime is non-default
    
    if flow_regime == FLOW_REGIME_LAMINAR:
        
        # laminar flow
        
        if (reynolds_number > LAMINAR_UPPER_RE_LIMIT or 
            reynolds_number < LAMINAR_LOWER_RE_LIMIT):
            
            # issue warning: re/flow regime mismatch
            
            Warning('The Reynolds number does not match the flow regime.')
            
            print('The Reynolds number does not match the flow regime.')
            
            # return with error
            
            return np.nan, np.nan
        
        elif reynolds_number == None:
            
            # midrange reynolds number
            
            reynolds_number = (LAMINAR_LOWER_RE_LIMIT+LAMINAR_UPPER_RE_LIMIT)/2
            
        #**********************************************************************
        #**********************************************************************
            
        # compute friction factor for laminar flow
            
        friction_factor = 64/reynolds_number
    
    else: # non laminar flows
    
        #**********************************************************************
    
        # the models are only valid within specified (re, rpr) limits

        # get the relative pipe roughness range for this model

        (model_rpr_lower, model_rpr_upper) = DFF_RPR_LIMITS[model]
    
        # get the reynolds number range for this model
            
        (model_re_lower, model_re_upper) = DFF_RE_LIMITS[model]

        # check if a reynolds number has been provided
        
        if reynolds_number == None:
            
            # no reynolds number has been provided: assume intermediate one
            
            reynolds_number = (model_re_lower+model_re_upper)/2
    
        if ignore_model_input_ranges == False:
            
            # check if the relative pipe roughness is within limits for this model
            
            if (pipe.e_eff/pipe.d_int > model_rpr_upper or
                pipe.e_eff/pipe.d_int < model_rpr_lower):
                
                # the relative pipe roughness is not within limits
                
                Warning('The relative roughness is out of range for this model.')
                
                print('The relative roughness is out of range for this model.')
                
                return np.nan, np.nan
            
            # check if a reynolds number has been provided
            
            if (reynolds_number < model_re_lower or 
                reynolds_number > model_re_upper):
                
                # reynolds number is not compatible with model range
                
                Warning('The Reynolds number is out of range for this model.')
                
                print('The Reynolds number is out of range for this model.')
                
                # return error
                
                return np.nan, np.nan
        
        #**********************************************************************
        #**********************************************************************
        
        # select how to proceed based on the model indicated
        
        if model == DFF_HAALAND:
            
            # Haaland, 1985
            
            friction_factor = DarcyFrictionFactorExplicitHaaland(
                reynolds_number, 
                pipe.e_eff, 
                pipe.d_int)
            
        elif model == DFF_PETHUKOV:
            
            # Petukhov, 1970
            
            friction_factor = DarcyFrictionFactorExplicitPethukov(
                reynolds_number)
            
        elif model == DFF_NIKURADSE:
            
            # Nikuradse, 1930
                
            friction_factor = DarcyFrictionFactorExplicitNikuradse(
                pipe.d_int,
                pipe.e_eff)
            
            # the reynolds number should be as high as possible for maximum
            # accuracy: 
                
            # override midrange point
            
            reynolds_number = MinimumReynoldsNumberCompleteTurbulence(
                friction_factor)
            
            # if reynolds_number > model_re_upper:
                
            #     reynolds_number = model_re_upper
            
        elif model == DFF_BLASIUS:
            
            friction_factor = DarcyFrictionFactorExplicitBlasius(
                reynolds_number)
            
        else: # unknown model
        
            Warning('Unknown Darcy friction factor model selected.')
            
            print('Unknown Darcy friction factor model selected.')
            
            return np.nan, np.nan
        
        #**********************************************************************
        #**********************************************************************
        
    return reynolds_number, friction_factor

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def LogMeanTemperatureDifference(temperature_inlet, 
                                 temperature_outlet, 
                                 temperature_surface):
    
    return (temperature_inlet-temperature_outlet)/(
        np.log((temperature_surface-temperature_outlet)/
               (temperature_surface-temperature_inlet))
        )

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def calculateFlowRegimeTuple(pipe: Pipe, 
                             fluid: Fluid,
                             fluid_speed=None, 
                             mass_flow_rate=None, 
                             volumetric_flow_rate=None):
    """
    Returns the fluid speed, the mass flow rate and the volumetric flow rate.

    Parameters
    ----------
    pipe : Pipe
        DESCRIPTION.
    fluid : Fluid
        DESCRIPTION.
    fluid_speed : TYPE, optional
        DESCRIPTION. The default is None.
    mass_flow_rate : TYPE, optional
        DESCRIPTION. The default is None.
    volumetric_flow_rate : TYPE, optional
        DESCRIPTION. The default is None.

    Returns
    -------
    fluid_speed : TYPE
        DESCRIPTION.
    mass_flow_rate : TYPE
        DESCRIPTION.
    volumetric_flow_rate : TYPE
        DESCRIPTION.

    """
    
    # at least one of the keyword arguments has to be non-default
    
    if (fluid_speed != None and
        mass_flow_rate == None and
        volumetric_flow_rate == None):
        
        # the fluid speed is known
        
        # volumetric flow rate
        
        volumetric_flow_rate = core.VolumetricFlowRateFromMeanFluidSpeed(
            fluid_speed, 
            pipe.InternalCrossSectionalArea())
        
        # mass flow rate
        
        mass_flow_rate = core.MassFlowRateFromVolumetricFlowRate(
            volumetric_flow_rate, 
            fluid.mass_density)
    
        #**********************************************************************
        #**********************************************************************
        
    elif (fluid_speed == None and
        mass_flow_rate != None and
        volumetric_flow_rate == None):
        
        # the mass flow rate is known
        
        # volumetric flow rate
        
        volumetric_flow_rate = core.VolumetricFlowRateFromMassFlowRate(
            mass_flow_rate, 
            fluid.mass_density)
        
        # fluid speed
        
        fluid_speed = core.MeanFluidSpeedViaVolumetricFlowRate(
            volumetric_flow_rate,
            pipe.InternalCrossSectionalArea())
    
        #**********************************************************************
        #**********************************************************************
    
    elif (fluid_speed == None and
        mass_flow_rate == None and
        volumetric_flow_rate != None):
        
        # the volumetric flow rate is known
    
        # fluid speed
        
        fluid_speed = core.MeanFluidSpeedViaVolumetricFlowRate(
            volumetric_flow_rate,
            pipe.InternalCrossSectionalArea())
        
        # mass flow rate
        
        mass_flow_rate = core.MassFlowRateFromVolumetricFlowRate(
            volumetric_flow_rate, 
            fluid.mass_density)
    
        #**********************************************************************
        #**********************************************************************
    
    elif (fluid_speed != None and
        mass_flow_rate != None and
        volumetric_flow_rate != None):
        
        # nothing to do here
        
        return fluid_speed, mass_flow_rate, volumetric_flow_rate
    
        #**********************************************************************
        #**********************************************************************
    
    else: 
        
        # there is a strange combination of inputs
        
        return np.nan, np.nan, np.nan
    
    #**************************************************************************
    #**************************************************************************
        
    return fluid_speed, mass_flow_rate, volumetric_flow_rate

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
    
def check_nusselt_number_model_is_valid(nusselt_number_model,
                                        reynolds_number: float, 
                                        prandtl_number: float,
                                        print_warnings: bool = False) -> bool:
    
    if nusselt_number_model not in NU_MODELS:
        
        raise NotImplementedError
        
    out = True
    
    if (reynolds_number < NU_RE_LIMITS[nusselt_number_model][0] or
        reynolds_number > NU_RE_LIMITS[nusselt_number_model][1]):
        
        if print_warnings:
        
            print('The Reynolds number (= ' + str(reynolds_number)+
                  ') is outside the valid range for the '+
                  str(nusselt_number_model)+
                  ' Nusselt number correlation. The valid range is between '+
                  str(NU_RE_LIMITS[nusselt_number_model][0])+
                  ' and ' +str(NU_RE_LIMITS[nusselt_number_model][1]) + '.')
        
        out = False
              
    if (prandtl_number < NU_PRANDTL_LIMITS[nusselt_number_model][0] or
        prandtl_number > NU_PRANDTL_LIMITS[nusselt_number_model][1]):
        
        if print_warnings:
        
            print('The Prandtl number (= ' + str(prandtl_number)+
                  ') is outside the valid range for the '+
                  str(nusselt_number_model)+
                  ' Nusselt number correlation. The valid range is between '+
                  str(NU_PRANDTL_LIMITS[nusselt_number_model][0])+
                  ' and ' +str(NU_PRANDTL_LIMITS[nusselt_number_model][1]) +
                  '.')
        
        out = False
          
    return out

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def NusseltNumber(pipe: Pipe,
                  friction_factor,
                  reynolds_number,
                  prandtl_number,
                  flow_regime,
                  laminar_condition = CONDITION_ISOTHERMAL_SURFACE,
                  turbulent_region_model = NU_GNIELINSKI,
                  ignore_reynolds_number_bounds: bool = True,
                  ignore_prandtl_number_bounds: bool = True):
    """Get the Nusselt number for internal forced convection in pipes."""
    
    if flow_regime == FLOW_REGIME_LAMINAR:
        
        # laminar flow
            
        if laminar_condition == CONDITION_CONSTANT_HEAT_FLOW:
            
            # constant heat flow condition
            
            return NusseltNumberLaminarEntranceConstantHeatFlow(
                pipe,
                reynolds_number,
                prandtl_number)
        
        elif laminar_condition == CONDITION_ISOTHERMAL_SURFACE:

            # isothermal condition
            
            return NusseltNumberLaminarEntranceIsothermal(
                pipe,
                reynolds_number,
                prandtl_number)
        
        else:
            
            # unknown laminar flow condition
            
            return math.nan
    
    elif flow_regime == FLOW_REGIME_TURBULENT:
        
        # turbulent
        
        # Cengel2003, page 443, says entrance regions for turbulent flow are
        # not particularly relevant
        
        if turbulent_region_model == NU_GNIELINSKI:
        
            check_nusselt_number_model_is_valid(turbulent_region_model,
                                                reynolds_number,
                                                prandtl_number,
                                                True)
        
            return _nusselt_number_gnielinski(
                reynolds_number, 
                prandtl_number,
                friction_factor)
        
        else:
            
            raise NotImplementedError
    
    else:
        
        # transition region
        
        # assume fluids behave as in the turbulence region
        
        check_nusselt_number_model_is_valid(NU_GNIELINSKI,
                                            reynolds_number,
                                            prandtl_number,
                                            True)
        
        return _nusselt_number_gnielinski(
            reynolds_number, 
            prandtl_number,
            friction_factor)

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def NusseltNumberLaminarEntranceIsothermal(pipe: Pipe,
                                           reynolds_number,
                                           prandtl_number):
    """Get the Nusselt number for laminar flow under the isothermal condition."""
    
    # Edwards et al., 1979
    
    return (3.66+
            0.065*(pipe.d_int/pipe.length)*reynolds_number*prandtl_number/
            (1+0.04*((pipe.d_int/pipe.length)*
                      reynolds_number*prandtl_number)**(2/3))
            )

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def NusseltNumberLaminarEntranceConstantHeatFlow(pipe: Pipe,
                                                 reynolds_number,
                                                 prandtl_number):
    """Get the Nusselt number for laminar flow under the constant heat flow condition."""
    
    return (4.36+
            0.065*(pipe.d_int/pipe.length)*reynolds_number*prandtl_number/
            (1+0.04*((pipe.d_int/pipe.length)*
                      reynolds_number*prandtl_number)**(2/3))
            )

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def _nusselt_number_gnielinski(reynolds_number, 
                               prandtl_number,
                               friction_factor):
    
    # Source: Cengel2003, page 441
    
    # Gnielinski (1976), Ref. 8
    
    # Primary source:    
    # V. Gnielinski. “New Equations for Heat and Mass
    # Transfer in Turbulent Pipe and Channel Flow.”
    # International Chemical Engineering 16 (1976), pp. 359–368.
        
    return (
        (friction_factor/8)*(reynolds_number-1e3)*prandtl_number/
        (1+12.7*math.sqrt(friction_factor/8)*(prandtl_number**(2/3)-1))
        )

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def NusseltNumber_ColburnDittusBoelter(reynolds_number, 
                                       prandtl_number,
                                       mode='colburn'):
    """Returns the Nusselt number using the Colburn or Dittus-Boelter equations."""
    
    # fully developed flow, turbulent        
    
    modes = ['colburn','heating','cooling']
    
    # Colburn/Dittus-Boelter equation according to Cengel (2003), page 441
    
    # Reynolds number:  > 10e3
    # Prandtl number: 0.7 <= Pr <= 160
    # Fluid properties should be evaluated at bulk mean fluid temperatures
    
    # Source: F. W. Dittus and L. M. K. Boelter. University of
    # California Publications on Engineering 2 (1930), p. 433.
    
    if mode not in modes:
        
        # raise UnknownOption
        
        return math.nan    
    
    if reynolds_number <= 10e3:
        
        # raise ReynoldsNumberOutOfRangeException
        
        return math.nan
        
    if prandtl_number < 0.7 or prandtl_number > 160:
        
        # raise PrandtlNumberOutOfRangeException
        
        return math.nan
    
    # get the nusselt number
    
    if mode == 'heating': # dittus-boelter equation for heating
        
        return 0.023*(reynolds_number**0.8)*(prandtl_number**0.4)
    
    elif mode == 'cooling': # dittus-boelter equation for cooling
        
        return 0.023*(reynolds_number**0.8)*(prandtl_number**0.3)
    
    else:
        
        # colburn equation
    
        return 0.023*(reynolds_number**0.8)*(prandtl_number**(1/3))

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def NusseltNumberRoughPipes(friction_factor,
                            reynolds_number,
                            prandtl_number):
    """Get the Nusselt number for internal forced convection in rough pipes."""
    
    # Cengel2003, page 441, equation 8-70
    # pr: 0.5 to 2000
    # re: 3e3 to 5e6
    
    if prandtl_number < 0.5 or prandtl_number > 2000:
        
        Warning('Prandtl number limits have not been observed.')
        
    if reynolds_number < 3e3 or reynolds_number > 5e6:
        
        Warning('Reynolds number limits have not been observed.')
        
    # return statement
        
    return ((prandtl_number*friction_factor/8)*(reynolds_number-1000)/
            (1+12.7*np.sqrt(friction_factor/8)*(prandtl_number**(2/3)-1)))

    # return (0.125*friction_factor*reynolds_number*prandtl_number**(1/3))
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************